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We present a model for the ion distribution near a charged surface, based on the response of the 
ions to the presence of a single test particle. Near an infinite planar surface this model produces 
the exact density profile in the limits of weak and strong coupling, which correspond to zero and 
infinite values of the dimensionless coupling parameter. At intermediate values of the coupling 
parameter our approach leads to approximate density profiles that agree qualitatively with Monte- 
Carlo simulation. For large values of the coupling parameter our model predicts a crossover from 
exponential to algebraic decay at large distance from the charged plate. Based on the test charge 
approach we argue that the exact density profile is described, in this regime, by a modified mean 
field equation, which takes into account the interaction of an ion with the ions close to the charged 
plate. 



I. INTRODUCTION 



dimensionless coupling parameter [2(j [31j 



Interactions between charged objects in solution are 
determined by the distribution of ions around them. Un- 
derstanding these distributions is thus of fundamental 
importance for theoretical treatment of water soluble 
macromolecules such as polyelectrolytes, charged mem- 
branes, and colloids 0,0- In recent years, much interest 
has been devoted to correlation effects in ionic solutions 
and to attempts to go beyond mean field theory in their 
treatment. In particular it has been realized that such 
effects can lead to attractive interactions between simi- 
larly charged objects, as was demonstrated in theoretical 
models |^MIM Irl [^lM[ ^.pjll, sim ulation [M lnllr3 . ll3LIT^ | 
and experiment [lllli lH 111 IS- 

Despite the theoretical interest in ion correlation ef- 
fects, they are not fully understood even in the most 
simple model for a charged object surrounded by its coun- 
terions, shown schematically in Fig. 1. The charged ob- 
ject in this model is an infinite planar surface localized 
at the plane z — 0, having a uniform charge density a. 
The charged plate is immersed in a solution of dielectric 
contact e and is neutralized by a single species of mobile 
counterions (there is no salt in the solution) . These coun- 
terions are confined to the half space z > and each one 
of them carries a charge e, which is a multiple of the unit 
charge for multivalent ions. The ions are considered as 
point-like, i.e., the only interactions in the system, apart 
from the excluded volume at z < 0, are electrostatic. 

The system described above is characterized by a single 
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where ksT is the thermal energy. At small values of 
this coupling parameter the electrostatic interaction be- 
tween ions is weak. As a result, in the limit 5 — > mean 
field theory is exact, as can be formally proved using a 
field-theory formulation of the problem [2£j . Correlations 
between ions close to the charged plate play a progres- 
sively more important role with increase of the coupling 
parameter. From Eq. Q one sees that this happens with 
an increase of the surface charge, with decrease of the 
temperature or dielectric constant, and with increase of 
the charge or, equivalently, the valency of counterions. 
The model of Fig. 1 thus provides an elementary theo- 
retical framework for studying ion correlation effects near 
charged objects, with no free parameters other than S, 
which tunes and controls the importance of ion correla- 
tions. 

In recent years two theoretical approaches were pro- 
posed for treatment of the strong coupling limit 5 — > 
oo. The first approach is based on properties of 
the strongly coupled, two dimensional one component 
plasma, and emphasizes the possibility of Wigner crystal 
like ordering parallel to the charged plane. The second 
approach 10] is formally an exact, virial type expansion 
applied to a field-theory formulation of the partition func- 
tion. Both of these approaches predict an exponential 
decay of the ion density distribution in the strong cou- 
pling limit, leading to a more compact counterion layer 
than in mean field theory. 

Although the form of the density profile is established 
in the two limits S — ► and S — > oo, its behavior at 
intermediate values of the coupling parameter is still not 
clear. Liquid-state theory approaches such as the AHNC 
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FIG. 1: Schematic description of the double layer model con- 
sidered in this work. An infinite charged plate having a uni- 
form surface charged density a is immersed in a dielectric 
medium having dielectric constant e on both sides of the plate. 
The charge of the plate is neutralized by point-like counteri- 
ons, carrying each a charge e. These ions are confined to the 
positive z half space, where z — is the plane occupied by 
the plate. In thermal (fcsT) units the interaction between 
two ions is given by Ib/t, where r is the distance between the 
ions and Ib = e 2 / '(efcsT) is the Bjerrum length. 



approximation |2l| can probably be used in this regime, 
but in practice they were applied in the literature only 
to relatively small values of the coupling parameter, usu- 
ally also including ingredients other than those in the 
model of Fig. 1 - such as finite ion size, added salt, or 
an interaction between two charged planar surfaces. The 
infinite planar double layer with no added salt (Fig. 1) 
was recently studied using Monte-Carlo computer simu- 
lation [l^, providing detailed results on the counterion 
distribution in a wide range of coupling parameter values. 
These results validate the expected behavior in the weak 
and strong coupling limits. In addition they provide new 
data at intermediate values of the coupling parameter, to 
which theoretical approaches can be compared. 

We propose, in the present work, a new theoretical ap- 
proach for treating the distribution of counterions near 
the charged plate. This approach is based on an ap- 
proximate evaluation of the response of the ionic layer, 
to the presence of a single test particle. While an ex- 
act evaluation of this response would, in principle, allow 
the distribution of ions to be obtained exactly, we show 
that even its approximate calculation provides meaning- 
ful and useful results. In the limits of small and large 
S the exact density profile is recovered. At intermedi- 
ate values of the coupling parameter our approach agrees 
semi-quantitatively with all the currently available sim- 
ulation data. 

The outline of this work is as follows. In Sec. |n] we 
present the model and discuss why it produces the exact 
density profile in the weak and strong coupling limits. 
In Sec. IIIII we present numerical results for the density 



profile close to the charged plate, and compare them with 
simulation results of Ref. [l3j . Numerical results of our 
model, further away from the charged plate, where there 
is currently no data from simulation, are presented in 
Sec. IIVI and scaling results are obtained for the behavior 
of our model in this regime. Finally, in Sec.|3we discuss 
the relevance of our model's predictions, at small and 
large z, to the exact theory. Many of the technical details 
and derivations appear in the appendices at the end of 
this work. 



II. MODEL 
A. Scaling 

Consider the system shown in Fig. 1, with the param- 
eters a, e, e defined in the introduction. We will first 
express the free energy using the dimensionless coupling 
parameter S. In the canonical ensemble the partition 
function can be written as follows (zi is the z coordinate 
of the z-th ion): 
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where lg — e 2 /ekBT is the distance at which the 
Coulomb energy of two ions is equal to the thermal en- 
ergy ksT, and fi = e/(2TrlB<j) characterizes the bare 
interaction of an ion with the charged plane. These 
quantities, the only two independent length scales in the 
problem, are known as the Bjerrum length and Gouy- 
Chapman length, respectively. We rescale the coordi- 
nates by the Gouy-Chapman length: 
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yielding exp(— F/v) = (/i) 3Af exp(— Fjy), where 
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and the ratio 
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is the coupling parameter that was previously defined in 
Eq. Q. The requirement of charge neutrality is: N/A = 
<r/e, where A is the planar area. Hence the number of 
ions per rescaled unit area is equal to: 
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where A = A//J 2 . The local density of ions in the rescaled 
coordinates is equal to p(r) = p?p{r). Due to symmetry 
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this density depends only on z and it is convenient to de- 
fine a normalized, dimensionless, density per unit length: 



n{z) = 2ttIb^ 2 P — 2n^p 



having the property: 



dzn(z) = 1 



(7) 
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as seen from Eqs. @l and J7J). From here on we will 
omit the tilde symbol (~) from all quantities, in order 
to simplify the notations. In order to express physical 
quantities in the original, non-scaled units, the following 
substitutions can be used: 
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We will also omit the subscript N from the free energy 
i*jv, implying that N is determined by charge neutrality. 



B. Known results in the limits of small and large S 

We briefly review some known properties of the ion 
distribution in the limits of small and large S (for a more 
complete discussion, see Ref. ^(J ) . In the limit of 3 — > 
mean field theory becomes exact. The density profile is 
obtained from the Poisson-Boltzmann (PB) equation and 
decays algebraically, having the form 



n PB (z) 



{z + lf 



(11) 



Within the adsorbed layer ions form a three dimensional, 
weakly correlated gas: the electrostatic interaction be- 
tween neighboring ions is small compared to the ther- 
mal energy. This last statement can be verified by con- 
sidering the density of ions at contact with the plane, 
Ppb(0) = 1/(2x3) (see Eqs. Q and (TTJl). The typical 
distance between neighboring ions is thus of order S 1 / 3 . 
In the non-scaled units this distance is much larger than 
Ib, which validates the statement that ions are weakly 
correlated: S 1 / 3 ^/ = S~ 2 / 3 ^ 3> Ib- Note also that this 
typical distance is small compared to the width of the 
adsorbed layer (Gouy-Chapman length): S 1 / 3 /! <C fi. 

In the opposite, strong coupling (SC) limit of S ^ 1, 
the density profile decays exponentially, 



nsc{z) = exp(-z) 



(12) 



The width of the adsorbed layer is still of order /i in 
the non-scaled units, but is now small compared to Ib- 
Equation indicates that the average lateral distance 
between ions is then of order S 1 / 2 . This distance is large 
compared to the width of the ionic layer, S 1 / 2 /! ^> p. On 
the other hand it is small in units of the Bjerrum length: 
S 1 / 2 /! = "Er x l 2 lB < !fl. The ions form, roughly speak- 
ing, a two-dimensional sheet and are highly correlated 



within this adsorbed layer. The typical lateral separa- 
tion between ions, 3 1//2 , is an important length scale in 
the strong coupling limit, and will play an important role 
also in our approximated model. 

At sufficiently large values of 5 it has been conjectured 
(but not proved) that ions form a two-dimensional, trian- 
gular close-packed Wigner crystal parallel to the charged 
plate. Based on the melting temperature of a two dimen- 
sional, one component plasma, one can estimate that this 
transition occurs at 5 > 31, 000 [HEl- Furthermore, the 
ion-ion correlation function is expected to display short 
range order similar to that of the Wigner crystal even 
far below this transition threshold. The exponential de- 
cay of Eq. Ill2[l was predicted, based on these notions, in 
Ref. |9j . The same result can be obtained also in a formal 
virial expansion [Toj|. which is valid for large 3 but does 
not involve long range order parallel to the charged plate 
at any value of H. 

Finally we note two general properties of the density 
profile that are valid at any value of 5. First, the nor- 
malized contact density n(0) is always equal to unity - 
a consequence of the contact theorem [23| (see also Ap- 
pendix [Dj . Second, the characteristic width of the ad- 
sorbed layer is always of order unity in the rescaled units. 
These two properties restrict the form of the density dis- 
tribution quite severely and indeed the two profiles i|ll|) 
and (|12J) are similar to each other close to the charged 
plane. Far away from the plate, however, they are very 
different from each other: at z > 1 the probability to 
find an ion is exponentially small in the SC limit, while 
in the weak coupling limit it decays only algebraically 
and is thus much larger. Furthermore, in the weak cou- 
pling case, moments of the density, including the average 
distance of an ion from the plate, diverge. 

Although the form of the density profile is known in the 
limits of small and large 3, two important issues remain 
open. The first issue is the form of the density profile at 
intermediate values of 5. At coupling parameter values 
such as 10 or 100 perturbative expansions around the 
limits of small or large 3 0,H(j are of little use, because 
they tend to give meaningful results only at small values 
of their expansion parameter. Such intermediate values 
are common in experimental systems with multivalent 
ions, as demonstrated in Table 1. Second, even at very 
small or very large 3 the respective asymptotic form of 
n(z) may be valid within a limited range of z values. 
In particular, for large 5 it is natural to suppose that 
sufficiently far away from the charged plate the density 
profile crosses over from SC to PB behavior. Indeed, far 
away from the plate the ion density is small, resembling 
the situation near a weakly charged surface. The main 
objective of this work is to investigate these issues, both 
of which necessitate going beyond the formal limits of 
vanishing and infinite 3. 
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TABLE I: Characteristic values of a, fi and S for several rep- 
resentative macromolecules. Values of S are shown for two 
cases: monovalent counterions (1-e) and 4-valent ones (4-e). 
The Gouy-Chapman length fj, corresponds to monovalent ions. 
The cell membrane charge density is estimated assuming that 
10 % of the lipids in the membrane are charged. The surface 
charge of DNA is estimated from the linear charge density 
along the DNA contour, equal to 1/1.7 e/A, and assuming a 
radius of 10 A. For Mica full dissociation of charged groups 
is assumed. In all three cases the Bjerrum length is taken as 
Ib = 7 A, which corresponds to water at room temperature. 
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C. Test-charge mean field model 

Our model is based on the following observation: the 
normalized density n(z) is proportional to the partition 
function of a system where a single ion is fixed at the 
coordinate z: 



n(z) = —exp[-F(z) 
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where the coordinate of the fixed (JV-th) ion in Eq. 114(1 
is ZqZ. Equations ((13|) - ((15(l are exact and can be readily 
formulated also in the grand canonical ensemble. 

In the original coordinates F(zo) is the free energy of 
ions in the external potential: 
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exerted by the charged plane and fixed ion. Examination 
of Eq. 1 (14(1 shows that in the rescaled coordinates these 
are ions of charge in the external potential: 



z + 



zqZ\ 



(17) 



Our starting point for evaluating n(z) is the exact re- 
lation expressed by Ea. ((13() but we will use a mean field 
approximation in order to evaluate F(zq). In this ap- 
proximation the free energy is expressed as an extremum 
of the following functional of <p: 



F PB (z ) 
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where ip is the reduced electrostatic potential, 0(z) is the 
Heaviside function, and i^ se if is an infinite self energy 
which does not depend on z . The derivation of Eq. 1 (18(1 
is given in Appendix lAl 

The mean field equation for ip is found from the re- 
quirement SFpB/S(p(r) = 0: 



X9(z)c- V - —6(z)+Z5(r- z z) (19) 

Z7T 



This equation describes the mean field distribution of 
ions in the presence of a charged plane of uniform charge 
density — l/(27r) (second term in Eq. ((19(1 ) and a point 
charge S located at r = zqz (third term in Eq. ((19(0 . In 
cylindrical coordinates the solution ip can be written as 
a function only of the radial coordinate r and of z, due 
to the symmetry of rotation around the z axis. 

It is easy to show that at the extremum of -Fpe the 
overall charge of the system, including the charged sur- 
face, test charge and mobile counterions, is zero. The 
fugacity A has no effect on the extremal value of Fp#; 
changing its value only shifts ip(r) by a constant. 

Equations ((13(1 and 1 (15(1 . together with the mean-field 
approximation for F(zq) given by Eqs. 1(18(1 and ((19(1 con- 
stitute the approximation used in this work: 



where 



n{z) = — exp[-F PB (z)] 
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We will refer to this approximation as the test-charge 
mean field (TCMF) model. 



D. Limits of small and large 5 

As a first example we present, in Fig. 2, density pro- 
files obtained numerically from Eqs. I(2(J|) - ((21|) at 3 = 0.1 
(circles) and at H = 10000 (squares). The continuous 
lines are the theoretically predicted profiles at 2 — * 0, 
^pb(z) = l/( z + I) 2 ; an( A at S ^ oo, nsc(z) = exp(— z). 
The figure demonstrates that the weak coupling and 
strong coupling limits are reproduced correctly in our 
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FIG. 2: Density profiles, n(z), numerically calculated using 
the TCMF model of Eqs. CSD-GD, with H = °- 1 (squares) 
and S =10,000 (circles). The solid lines show the exact 
asymptotic profiles in the low coupling, jipb(z) = 1/(2 + l) 2 , 
and in the strong coupling limit, nsc = exp(— z). 



approximation. Before presenting further numerical re- 
sults, we discuss the behavior of our model in the two 
limits of small and large S. 

Our discussion is based on the following exact identity: 



dz 



log [n(z )] 



dz 



d 

F(z ) = - -- W>(r;«o)) 



= z z 
(22) 

where (ip(r;zo)) is the thermally averaged electrostatic 
potential at r, when a test charge is fixed at zqZ (the first 
argument of (<p(r; Zq)} designates the position r where the 
potential is evaluated, while the second argument des- 
ignates the position of the test charge, z z). In other 
words, the gradient of log[n(zo)] is equal to the electro- 
static force acting on a test charge positioned at r = zqz. 
This equation does not involve any approximations and 
is proved in Appendix iBl 

Within our approximation, where F(zq) is replaced by 
-Fpb(zo), a similar equation holds (also proved in the Ap- 
pendix) : 



-^logKzo)] = -^F PB (zo) = - 
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dz 



= z z 
(23) 

where (p(r\zo) is now the solution of Eq. 1(19(1 . In other 
words, the gradient of log[n(z )] is equal to the electro- 
static force experienced by a test charge positioned at 
r = z z, evaluated using the mean field equation (|19|l . 
This quantity, 



f(z ) 



dz 



(24) 
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will be studied in detail below because of its important 
role within our model. With this notation the relation 



between f(z) and n(z) reads: 
d 



dz 



log [n(z)] = -/(*) 



(25) 



Using Eq. i(25|) we can understand why both the weak 
and strong coupling limits are reproduced correctly in 
our model: 



Weak coupling 



In the limit a — ► 0, 



-^-ip(r;z ) -> ^-(pp B {z). 



(26) 



where tfp B (z) is the solution of Eq. 1)19(1 without a test 
charge, i.e., setting S = 0. We note that the potential (p 
(Ea. 1191) has three sources: the charge of mobile counteri- 
ons, X8(z)e~ ip , the uniformly charged plane, and the test 
charge. Although the potential due to the test charge is 
infinite at r = zqz, its derivative with respect to z is zero 
and has no contribution in Eq. 1(26(1 . Using Eq. 1(25(1 we 
find: 



-^-log [n(z)] = ~^-tp PB (z) 
dz dz 



(27) 



This equation, together with the normalization require- 
ment for n(z) leads to the result: 

n(z) = — exp[-<£ PB 0z)] = n PB {z) (28) 
^0 



Strong coupling 

In the strong coupling limit, 5 — ► 00, a correlation hole 
forms in the distribution of mobile counterions around 
the test charge at r = z$z. The structure and size of 
this hole, as obtained from Eq. ((19(1 . will be discussed 
in detail later. For now it is sufficient to note that 
the correlation hole gets bigger with increasing H. As 
5 — > 00 the force at z$z due to the mobile counterions 
vanishes, leaving only the contribution of the charged 
plane: (d/dz)tp(r; Zq, H) — * 1. Hence in this limit 



— log[n(2)] = -l 
dz 

leading to the strong coupling result: 
n(z) — exp(— z) 



(29) 



(30) 



where the prefactor of the exponent follows from the nor- 
malization condition, Eq. (jSJ). 

In the rest of this work we will explore predictions of 
the TCMF model at intermediate coupling, where nei- 
ther of the two limits presented above is valid. Before 
proceeding we note that a similar discussion as above, of 
the weak and strong coupling limits, applies also to the 
exact theory, because of Eq. I(22|) . 
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III. NUMERICAL RESULTS AND 
COMPARISON WITH SIMULATION 

Results for f(z) 

We consider first the behavior of f(z), defined in 
Eq. (Pljl . close to the charged plate. Fig. 3(a) shows 
this behavior for 5 = 1, 10, 10 2 , 10 3 and 10 4 (alternating 
solid and dashed lines). The curves were obtained from 
a numerical solution of the partial differential equation 
(PDE), Eq. O (see Appendix IU1 for details of the nu- 
merical scheme). For comparison the weak coupling (PB) 
and strong coupling (SC) limits are shown using dotted 
lines: 

/pb(z) = 7TT ; f*°& = 1 ( 31 ) 

As S increases f(z) gradually shifts from PB to SC be- 
havior. At E = 10 4 , f(z) is very close to 1 within the 
range of z shown in the plot, although there is still a 
noticeable small deviation from unity. 

In Fig. 3(b) these results are compared with simula- 
tion data (symbols), adapted from Ref. 01 . The value 
of f(z) was obtained from the simulation results for n(z) 
using the relation dlog[n(^)]/d2 = —f(z) 32]. Qualita- 
tively our results agree very well with simulation. Note 
especially the gradual decrease of f(z) with increasing z 
for S = 100 (diamonds): this value of S is far away from 
both the weak coupling and the strong coupling limits. 
The regime where f(z) decreases linearly with z will be 
further discussed in Sec. IIV Al 

It was previously conjectured [l(j that for all values of 
S the SC limit is valid close enough to the charged plane. 
We note, however, that at contact with the plane f(z) is 
different from unity at small and intermediate values of 
5. Hence it is not very meaningful to define a region 
close to the plane where the SC limit is valid, unless S 
is very large. Values of f(z), extracted from simulation 
data in Fig. 3(b), suggest the same conclusion, i.e., f(z) 
does not approach unity at contact with the plane. A 
more accurate measurement of f(z) in the simulation is 
desirable because the error bars, as obtained in Fig. 3(b), 
are quite large. 

Results for n(z) 

The density profile n(z) can be found numerically by 
integrating Eq. (|25|) and use of the normalization condi- 
tion (jSJl 33]. Figure 2 already demonstrated that n{z) 
coincides with the exact profiles, 7Jpb(z) and nsc( z )> m 
the limits of small and large 2. Figure 4(a) shows the 
difference between n(z) and upb{z) for S = 1, 10, 10 2 , 
and 10 4 , as calculated numerically in the TCMF model 
(continuous lines). These results are compared with sim- 
ulation data (symbols) (l3L |24| . 

We first observe that the contact theorem is not obeyed 
in our approximation: n(0) — ?i.pb(0) = n(0) — 1 is differ- 



ent from zero. This is an undesirable property, because 
the contact theorem is an exact relation. The contact 
theorem is obeyed in the TCMF model only in the lim- 
its of small and large S, where the density profile as a 
whole agrees with the exact form, and the normalization 
condition JSJ enforces n(0) to be correct. The violation 
at intermediate values of S is finite, small compared to 
unity, and peaks at 3 between 10 and 100. At these val- 
ues of S the overall correction to PB is quite inaccurate 
at the immediate vicinity of the charged plate. On the 
other hand, at z larger than 1 our approximated results 
agree quite well with simulation for all the values of 5, 
as seen in Fig. 4(a). 

The violation of the contact theorem in the TCMF 
model can be traced to a non-zero net force exerted by 
the ionic solution on itself (see Appendix IP)l . This incon- 
sistency results from the use of a mean field approxima- 
tion for the distribution of ions around the test charge, 
while the probability distribution of the test charge itself 
is given by Eq. (|2Ti|) . 

It is possible to evaluate exactly the first order correc- 
tion in S to the PB density profile in the TCMF model, 
the details of which are given in Appendix El This cor- 
rection turns out to be different from the exact first order 
correction, which was calculated in Ref. |20J using a loop 
expansion up to one loop order (also shown in the Ap- 
pendix). It is important to note, however, that the lat- 
ter correction provides a useful result only for relatively 
small values of S. At S of order 10 and larger TCMF 
results are much closer to simulation than those of the 
loop expansion. 

Further comparison with simulation is shown in 
Fig. 4(b). Here we show the density n(z) itself, rather 
than the difference with respect to npsiz). The data 
is shown on a larger range of z than in part (a) and a 
logarithmic scale is used in order to allow small values of 
n(z) to be observed far away from the plate. In the range 
shown the TCMF model agrees semi-quantitatively with 
simulation. 

As a summary of this section we can say that the test 
charge mean field (TCMF) model captures the essential 
behavior of the ion distribution at close and moderate 
distances from the charged plate. Furthermore, all the 
available data from simulation agrees qualitatively with 
our approximation's predictions. 



IV. TCMF RESULTS FAR AWAY FROM THE 
CHARGED PLATE 

Our analysis of the ion distribution far away from the 
charged plate is done, at first, strictly within the context 
of the TCMF model, while a discussion of its relevance to 
the exact theory is deferred to Sec.EJ The main question 
of interest is whether a transition to PB behavior occurs 
at sufficiently large z, even for large values of 2. 

As a first step we will identify the important length 
scales characterizing the density distribution. Let us con- 
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Z Z 

FIG. 3: (a) Electrostatic force acting on a test charge, f(z), numerically calculated using the mean field equation 119^ . 
Alternating solid and dashed lines show f(z) for H = 1, fO, f0 2 ,f0 3 , and fO 4 . The dotted lines show the PB and SC forms 
of f(z), fpB(z) = l/(z + 1) and fsc(z) = exp(— z). (b) Comparison of f(z), calculated from Eq. (| 1 9|l (solid and dashed lines, 
same as in part (a)), with results from Monte-Carlo simulation |2^|. adapted from Ref. [l3| (H = 1, circles; H = 10, squares; 
5 = 10 2 , diamonds; 3 = 10 4 , triangles). Values of f(z) are obtained from simulation results for n(z) using the exact relation 
dlogn(,z)/dz = —f(z). Numerical estimation of the derivative of log [n(z)] results in relatively large error bars, which are shown 
as vertical lines. 



centrate first on the size of the correlation hole around a 
test charge. Naively we may expect this size to be of or- 
der 2, due to the form of the bare potential, S/|r — z n z\. 
A simple argument shows, however, that when the test 
charge is close to the charged plane the size of the cor- 
relation hole is much smaller than EE. Assume, roughly, 
that the mobile ion density is zero within a cylindrical 
shell of radius R around the test charge. The amount of 
charge depleted from this cylinder is then equal to i? 2 /2, 
since the surface charge per unit area is equal to 1/27T 
(see Eq. flT9)). This depleted charge must balance ex- 
actly the charge of the test particle, equal to EE, yielding 
a cylinder radius that scales as \/S rather than as EE. 

In order to put this argument to test, Fig. 5(a) shows 
the density of mobile ions calculated from Eq. I|19fl with 
a test charge at zq = 1, having H = 1,000. The shape 
of the correlation hole is roughly cylindrical and its ra- 
dius is indeed of order Vs ~ 30. The influence of the 
test charge on its surroundings is very non-linear, with a 
sharp spatial transition from the region close to the test 
charge, where the density is nearly zero, to the region 
further away, where the effect of the test charge is very 
small. At larger separations from the plate the quali- 
tative picture remains the same, as long as zq is small 
compared to vEE and provided that E>1. 

A very different distribution of mobile ions is found 
when zo is of order EE, as shown in Fig. 5(b). The cou- 
pling parameter is the same as in part (a), E = 1,000, 
but the test charge is now at z — 1,000. Instead of 
showing directly the density of mobile ions as in part 
(a), the figure shows the ratio between this density and 
^pb(z) = l/(z + l) 2 . This ratio is now very close to 
unity near the charged plane, where most of the ions 



are present. It is small compared to unity only within a 
spherical correlation hole around the test charge, whose 
size is of order EE. 

The above examples lead us to divide our discussion of 
the z dependence into two regimes: 

A. z < VE 

In order to justify use of the cylindrical correlation 
hole approximation within this range, let us assume such 
a correlation hole and calculate the force acting on the 
test charge: 

/(*)- r°d^ PB (*')— z '~ z (32) 

Jo V E +( Z - z ) 

where 71pb(z') is given by Eq. pi|l . the radius of the 
cylindrical region from which ions are depleted is taken 
as \/S, and the expression multiplying tipb(z') is the 
force exerted by a charged sheet having a circular hole 
of radius y/E, and positioned in the plane z'. Figures 
6 (a) and (b) show a comparison of this approximation 
(dashed lines) with that obtained from a full numerical 
solution of Eq. I|19|) (solid lines). The coupling parameter 
is equal to 10,000 in (a) and to 100 in (b). In both cases 
the approximation works well up to zq ~ In Fig. 6 
(c) the force acting on a test charge at contact with the 
plane, /(0), is shown for five values of EE! between 1 and 
10,000 (symbols), and compared with the approximation 
of Eq. J3IjJ (solid line). Note that Eq. is not a good 
approximation when EE is of order unity or smaller, since 
the correlation hole is then small compared to the width 
of the ion layer. 



8 



Q - - 


- H=1 


■-■ 


- 3=10 




3=100 


□ — 


S=1 0,000 




Ref. M, 



10" 



10" 
c' 10" 




10" 4 




■-■ 


- S-10 






3=100 


10" 5 




A 


3=1000 



10 



10 



FIG. 4: (a) Correction to the PB density profile, n(z) — 
npB(z), calculated numerically using the TCMF model, as 
function of z (lines). For comparison, symbols show the cor- 
rection obtained from Monte-Carlo simulation ^j, 24]. Four 
values of H are shown (see legend), 1, 10, 10 2 , and 10 4 . (b) 
The density profile itself, n(z), on a wider range of z than 
in part (a) and using logarithmic scales in both axes (lines - 
TCMF; symbols - MC simulation). 



B. 



z > Va 



When the test charge is far away from the plane, its 
effect close to the charged plate is small, suggesting that 
a linear response calculation may be applicable: 



f{z)=fp B {z)+Ef l (z) 



(33) 



The first term in this equation is the PB value of f(z), 
while fi(z) can be calculated using previous results of 



1 dg(z) 



2 dz A{z + l) 3 
- (l+ije^ 1 -*)* [l-z + (l-2i)z 2 + z 3 ] 

xE t [(l-i)z] 
- (1 - i)e (1+t)z [1 - z + (1 + 2i)z 2 + z 3 } 

x£i[(l +*)*]} (34) 

where g(z) was defined in Ref. pfij and is given by 
Eq. I)E11|) . and E\{z) is the exponential- integral func- 
tion [25j . We prove the first equality of Eq. I|34|) in Ap- 
pendix IE1 Figure 7 shows fi{z) (solid line) together with 
its asymptotic form for large z (dashed line), 



z > 1 



(35) 



Note that the asymptotic form of S/i(z) has the same z 
dependence as the electrostatic force exerted by a metal- 
lic surface, equal to S/(4z 2 ) in our notation, but the 
numerical prefactor is different. 

Although the influence of the test charge is small near 
the surface, its influence on ions in its close vicinity is 
highly non-linear and definitely not small. Hence the ap- 
plicability of Eq. (l-i-il) is far from being obvious when 
3 is large. We check this numerically by calculating 
f(z) — /pb( 2: ) an d comparing with S/x(z). The results 
are shown in Fig. 8(a), for 5 values of 3: 1, 10, 10 2 , 10 3 , 
and 10 4 . 

For each value of 3 the ratio (/ — /pb)/(3/i) (shown 
in the plot) approaches unity as z is increased, showing 
that Eq. I)33|l does become valid at sufficiently large z. 
The approach is, however, rather slow: a value close to 
unity is reached only when z>S. At z — 3 the ratio is 
approximately equal to 0.6 in all five cases. We conclude 
that the linear approximation of Eq. (|33|l is applicable 
only for z » S. Note that at these distances from the 
charged plate the linear correction itself is very small 
compared to the PB term, 



5A(z) 35 z + l 
f PB (z) ~ 4z 2 2 



If*' 



(36) 



where we also assumed that z 3> 1 and used Eq. (|35|l . 

Further insight on the results shown in Fig. 8(a) is ob- 
tained by noting that all of them approximately collapse 
on a single curve after scaling the z coordinate by 3. This 
curve, denoted by /i(z/5), is shown in Fig. 8(b): 



f(z)-f PB (z)^Ef 1 (z)xh^) 



(37) 



In order to demonstrate at what range of z values this 
scaling result is valid the same data is shown in Fig. 9 
using a logarithmic scale in the horizontal (z/5) axis. It 
is then seen clearly that (|37|l holds for z/3 larger than 
a minimal value, which is proportional to S -1 / 2 . The 




FIG. 5: (a) Scaled density of ions around a test charge, positioned at zq — 1, as obtained from Eq. Ijf 9|l . The cross designates 
the position of the test charge. In cylindrical coordinates the density is a function only of z (horizontal axis) and r (vertical 
axis). Darker shading in the plot means larger density (see also the legend on the right). The coupling constant is 3 = 1,000. 
For r larger than \/3 the profile, as function of z, quickly converges to the PB profile, upb(z) (b) A similar plot as in part (a), 
but the test charge is now at zq = 3 = 1,000. Here the ratio between the density and the PB profile is shown, rather than the 
density itself. This ratio is everywhere a number between zero (black) and one (white). The effect of the test charge on the 
ion distribution is large only within a correlation hole around the test charge, having approximately a spherical shape and a 
radius of order 3. 



vertical arrows designate z/S = 1.5/%/S for each value 
of S, approximately the point where the scaling becomes 
valid. Returning to consider z itself, we conclude that 
(|37|l holds for z > 1.5-\/S. This result justifies the sep- 
aration of the z dependence into two regimes, z < \/S 
and z > vH- 

We finally turn to consider the ion density n(z). Using 
Eqs. and lj»T7|) we can write 

-A n(z) = f {z) ~ /pB ( z ) + |-| (|)" 2 ^ (I) (38) 

which leads, upon integration, to the approximate scaling 
result, 

n{z) ~ c (S) * r, (|) (39) 

where 

^[r ) {u)} = -\ f°°du'^ (40) 
4 J u u' 2 

The prefactor cq(5) depends, through the normalization 
condition, on the behavior of n(z) close to the charged 
plane where the scaling form of Eq. (|39|l is not valid. 
Equation l|39(l is indeed validated by the numerical data 
and applies for z > \/S and S ^> 1. 

The density itself, n(z), is plotted in Fig. 10 using a 
semi-logarithmic scale, for S = 1, 10, 10 2 , 10 3 , and 10 4 . 
At z > S n(z) is proportional to ripe(z), as expected 
from Eq. (|39fl . The prefactor Co is extremely small for 
large S. We recall that Co is mainly determined by the 
behavior close to the charged plane, where f(z) is of order 
unity up to z < Following this observation we can 



expect log[co(H)] to scale roughly as This estimate is 
validated by the numerical results and is demonstrated 
in the figure using the horizontal arrows. For 5 = 10, 
10 2 , 10 3 , and 10 4 these arrows show an estimate for n(z) 
at z — 10 4 , given by 

upb(z) x exp(— 0.8VS) 
in very good agreement with the actual value of n. 



V. FURTHER DISCUSSION 

At this point we may ask to what extent our results 
represent the behavior of n(z) in the exact theory. Be- 
fore discussing this question we turn our attention for 
a moment to the system that our problem was mapped 
into in Eq. H19|) - that of a single ion of valence S in 
contact with a charged plane and its monovalent countc- 
rions. The results of Sections ITTT1 and Hvl can be regarded 
as exact for such a system, provided that the monovalent 
ions are weakly correlated (having, by themselves, a small 
coupling parameter as determined from their charge and 
that of the planar surface) . These results are thus of di- 
rect relevance to the interaction of a large multivalent 
macroion with a charged surface that is immersed in a 
weakly correlated solution of countcrions. 

Returning to the original question, we separate our dis- 
cussion according to the scaling results of the numerical 
analysis: 
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FIG. 6: (a), (b) Comparison of the approximation to f(z) 
given by Eq. 1321 (dashed line) , with a full numerical solution 
of the PDE (solid line). The coupling parameter E is equal to 
10,000 in (a) and to 100 in (b). Note that the approximation 
shown in the dashed line is good up to a distance from the 
plate equal to about \/3 in both cases. A distance of \/H from 
the charged plate is designated by the vertical dotted lines. 
Part (c) shows a comparison of /(0) in the approximation 
given by Eq. 1321 (symbols) and in the exact PDE solution 
(solid line) for a wide range of H values. 



z < V3 



N 




FIG. 7: First order (linear) term in an expansion of f{z): 
f(z) — /pb(z) + H/i(z) + Eq. 11341 . as obtained from 
the loop expansion of Ref. [2(J. The dashed line shows the 
asymptotic form of fi{z) at large z, fi(z) — 3/(4a 2 ). 



is obtained in our model, as well as in simulation and in 
the perturbative strong coupling expansion of Ref. |10| . 
However, at intermediate values of S such as 10, 100, 
or 1000 our results show that this exponential decay is 
only a rough approximation. The decoupling of a test 
charge from the rest of the ions is only partial, even at 
z = 0, leading to a value of f(z) that is (i) larger than 
unity at z = and (ii) considerably smaller than unity 
at z = vH- Both of these predictions are validated by 
simulation, as shown in Fig. 3(b). 

The quantitative agreement in f(z) between our model 
and simulation is surprisingly good, considering that the 
ionic environment surrounding the test charge is differ- 
ent in our approximation, compared to the exact the- 
ory. This good agreement can be attributed to the cor- 
rect length scales that characterize the approximate ionic 
environment: VS in the lateral direction, and 1 in the 
transverse direction. Indeed, in the lateral direction our 
results can be compared with pair distributions that were 
obtained in Monte-Carlo simulation The pair dis- 

tributions found in the simulation are characterized by a 
strong depletion within a correlation hole having diame- 
ter of order VS, in great similarity to Fig. 5(a). What 
is not captured by our approach is that multiple maxima 
and minima exist at S > 100 j3j|. Nevertheless, even at 
the very large coupling parameter value E = 10 4 these 
oscillations decay quite rapidly with lateral distance, and 
we can still say that the TCMF model captures the most 
significant feature of the pair distribution (namely, the 
structure of the correlation hole). 



When 5 is very large a test charge at z < is es- 
sentially decoupled from the rest of the ionic solution, 
feeling only the force exerted by the charged plane. This 
is the reason why an exponential decay, n[z) ~ exp(— z), 



B. 



z > 



Throughout most of this section we concentrate on the 
case z > 5, while a short discussion of the range < 
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FIG. 8: Comparison between the correction to f(z) relative to /pb(z), with the linearized expression E/i(z). In (a) the ratio 
[f(z) — /pb(z)]/[H/i(z)] is shown as function of z for five different values of H: 1, 10, 10 2 , 10 3 , and 10 4 (see legend in part (b); 
symbols show the same quantity as the solid line and are displayed in order to distinguish between the five lines). The ratio 
approaches unity at z much larger than H. In (b) the same data as in (a) is shown as function of z/H, leading to an almost 
perfect collapse of the five data sets on a single curve. 
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FIG. 9: Same data as in Fig. 9(b), shown using a logarithmic 
scale in the horizontal (z/H) axis. The approximated collapse 
of the different data sets, corresponding to different values of 
H, is seen to be valid only in the regime z > \/H. The vertical 
arrows mark z = 1.5yS for 3 = 1, 10, 10 2 , and 10 3 , above 
which the scaling @ is approximately valid. 



z < S is presented at the end of this section. 

Our model predicts a transition to algebraic decay of 
n(z) at z > S. Similar predictions were made in Rcf. 9] 
and in Ref. 0], where it was estimated that mean field 
results are valid for z > SlogS based on a perturbative 
expansion around mean field. There are currently no 
available results from simulation in this regime. 

A mean field behavior is obtained in our model in the 
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z 

FIG. 10: Scaled ion density, n(z) calculated using the TCMF 
model, shown for five different values of S (solid lines, top to 
bottom): 1, 10, 10 2 , 10 3 , and 10 4 . A logarithmic scale is 
used on both axis, allowing the behavior far away from the 
charged plate to be seen. The dashed lines show npe (z) (up- 
per dashed line) and nsc(z) (lower dashed line). At z > S 
the density profile is proportional to npa(z), with a prefactor 
whose logarithm scales as vH- To demonstrate this the hori- 
zontal arrows mark the value of exp(— 0.8\/H) x npp,(10 4 ) for 
S = 10, 10 2 , 10 3 , and 10 4 . The prefactor 0.8 is an approxi- 
mate fitting parameter. 

sense that 

/(z)^/ PB (*) = _L_ (41) 

decays as 1/z for large z. Nevertheless the finer details 
in our results do not match the form predicted by PB 
theory. The starting point for the following discussion 
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is an hypothesis that sufficiently far from the plate the 
exact density follows a PB form, 



^-asyrn 



1 



{z + bf 



(42) 



where b (or fib in the original, non-scaled coordinates) 
is an effective Gouy-Chapman length, characterizing the 
ionic solution far away from the plate. 

The asymptotic density profile found in our approxi- 
mation, n(z) — c(S)/(z + l) 2 , is different from Eq. 1|42JI . 
To understand this difference let us explain first the 
asymptotic behavior of f(z); it decays in the TCMF 
model as l/(z+ 1) because beyond the correlation hole 
the test charge is surrounded by an ion density of the 
form ripe = I/O? + l) 2 - This form is different from the 
profile n(z) that is eventually obtained by integrating 
f(z) - an inconsistency which is the source of the differ- 
ence between /(z) in our approximation and the hypoth- 
esized form f(z) ~ l/(z + b) (see also the discussion in 
Appendix [D]) . 

The behavior of our approximate f(z) leads to a decay 
of n(z) of the form l/(z + l) 2 - The normalization condi- 
tion for n(z) is then enforced through a small prefactor 
c(S). In comparison, in the hypothesized form (|42|l the 
prefactor must be 1 and the normalization is achieved by 
a large value of b. Note that b must be an extremely large 
number for large values of S: due to the exponential de- 
cay of n(z) at z < \/S the logarithm of b must be at least 
of order yE. 

Further insight on the behavior at z > S can be ob- 
tained using the exact equation (|22|l : 



^ log [n(z)] = -f(z) 



(43) 



where f(z) is now the mean (thermally averaged) elec- 
trostatic force acting on a test charge at distance z from 
the plate, in the exact theory. 

For the mean field form n a sym(z) to be correct, the 
contribution to f(z) coming from the influence of the test 
charge on its environment must be small compared to the 
mean field force, which is given by l/(z + b). Following 
our results from the previous section, the former quantity 
is given by qH/z 2 , where a is of order unity. Using the 
mean field equation l|19l) we obtained a = 3/4; in the 
exact theory, and for large H, where ions form a much 
more localized layer close to the plane than in mean field, 
it is plausible that a — 1/4, as in the force acting on a 
test charge next to a conducting surface 0, ^| . In any 
case, for Eq. (|42|l to represent correctly the decay of n(z) 
we must have 



oS 1 

—7T < 



z + b 



(44) 



leading to the result z > (S6) 1 / 2 which is exponentially 
large due to b. Up to this crossover distance the decay of 
n(z) is dominated by the correlation-induced interaction 
with the ions close to the plate. 



We conclude that for a very large range of z values 
the decay of n(z) must be different from Eq. I|42|l . At the 
same time, a mean field argument is probably applicable, 
because the density of ions is very small in this regime: 
we may presume that the contribution to /(z) can be 
divided into two parts - one part, coming from ions far 
away from the plate, which is hardly influenced by the 
test charge; and a second part, coming from ions close to 
the charged plate, where the test charge influence on f(z) 
is given by aS/z 2 . This reasoning leads to the following 
differential equation for n(z): 



d 2 , r . S1 , . 2aZ 
— log[n(z)] = 2n[z) + — 



(45) 



whose detailed derivation is given in Appendix [FJ By 
defining n(z) — exp(— (f> + aE/z) Eq. ll4*5|l is recast in the 
form: 



-^ = ~2exp[-0+^-=) =-2n(z) (4(>) 



showing that mean field theory is applicable, but an ex- 
ternal potential — otE/z, coming from the ions close to 
the plate, must be included in the PB equation. In prac- 
tice, for large S this equation will lead to a decay of the 
form n(z) ~ exp(aS/z) as suggested also in Refs. |9Hl3j. 
while a crossover to an algebraic decay will occur only 

I 1 / 2 ™ 
where Eq. l|44|) 



at a distance of at least 



Sexp(v / S) 

has been used and prefactors of order unity, inside and 
outside the exponential, are omitted. A numerical solu- 
tion of Eq. 14(j|) may be useful in order to describe the 
ionic layer at intermediate values of S (of order 10 - 100), 
where both mean field and correlation-induced forces are 
of importance at moderate distances from the plate. In 
order to test this idea quantitatively more data from sim- 
ulation is required. 

Finally we discuss the case where z > y/S but z is 
not large compared to S. Let us also assume that S 
is very large, so that most of the ions are much closer 
to the plate than a test charge fixed at zz. Within the 
TCMF model the effect of the test charge on ions close 
to the plate is nonlinear, leading to the scaling result of 
Eq. (|37|l . Similarly, in the exact theory it is not clear 
whether the correlation-induced force acting on the test 
charge is of the form aS/z 2 , since the effect of the test 
charge on ions close to the plate is not a small pertur- 
bation. Hence we believe that the relevance of Eq. I|46|) 
for z < 5, and the behavior of /(z) in this regime, merit 
further investigation [35|. 



VI. CONCLUSION 

In this work we showed how ion correlation effects can 
be studied by evaluating the response of the ionic solu- 
tion to the presence of a single test particle. Although 
we calculated this response using a mean field approxi- 
mation, we obtained exact results in the limits of small 
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and large 3, and qualitative agreement with simulation 
at intermediate values. 

The approach taken in this work demonstrates that 
for highly correlated ionic liquids it is essential to treat 
the particle charge in a non-perturbative manner. Once a 
single ion is singled out, even a mean field approximation 
applied to the other ions provides useful results. This 
scheme, called the test-charge mean field (TCMF) model, 
provides a relatively simple approximation, capturing the 
essential effects of strong correlations - to which more 
sophisticated treatments can be compared. 

Technically what is evaluated in this work is the ion- 
surface correlation function. Consideration of correlation 
functions of various orders leads naturally to liquid state 
theory approximations 26], some of which are very suc- 
cessful in describing ionic liquids [2l| . In particular these 
approximations usually treat the ion-ion correlation func- 
tion in a more consistent manner than the approximation 
used in this work, thus possibly alleviating some of the 
undesirable features of the TCMF model presented here, 
such as the violation of the contact theorem. The main 
advantage of the TCMF model is its simplicity, allowing 
the behavior of n(z) to be understood in all the range 
of coupling parameter values in terms of f(z), the force 
acting on a test charge. Furthermore, the exact equation 
Q22JI. which does not involve any approximation, is a use- 
ful tool in assessing correlation effects - as was done, for 
example, in this work in the end of Sec. 

It will be useful to summarize the main results ob- 
tained in this work. First, the exact equation 11' I'll pro- 
vides a simple explanation of the exponential decay of 
the density profile in the strong coupling limit. In light 
of this equation, exponential decay is expected as long 
as the test charge is decoupled from the rest of the ionic 
solution. Note that there is no necessity for long range or- 
der to exist in order for the exponential decay to occur, 
as was emphasized also in Ref. 9j. Indeed, within our 
TCMF model the ion distribution around a test charge 
does not display any long range order (see Fig. 5(a)) and 
yet simulation results, in particular in the strong coupling 
limit, are recovered very successfully. 

Second, the characteristic size of the correlation hole 
around an ion close to the plane, v3j plays an important 
role in determining the density profile. For very large 5 
the profile decays exponentially up to z < \/3, beyond 
which a crossover to a less rapid decay occurs. For inter- 
mediate values of the coupling parameter z = \/S is still 
an approximate boundary between regimes of different 
behavior of n(z), but the density profile at z < V3 does 
not decay in the simple exponential form exp(— z). In 
this sense one cannot speak of a region close to the plate 
where strong coupling results are valid. 

For z > 3 our approximation predicts a transition to 
an algebraic decay of n(z), of the form c(z)/(z + l) 2 , 
where the prefactor c(z) is exponentially small for large 
3. A different asymptotic behavior of the form l/(z + b) 2 
is probably valid at very large z, but is not predicted by 
the TCMF model. Arguments presented in Sec.0 based 



on the exact equation (1221) , lead to the conclusion that for 
large 3 the latter form (with a constant value of 6) can 
only be valid at extremely large values of z, while sug- 
gesting that at all distances from the plate larger than 2 
a modified mean field equation, Eq. (|H)|> . is valid. This 
equation, matched with the behavior of the ion distribu- 
tion close to the charged plate, ultimately determines the 
value of the effective Gouy-Chapman length b. 

Finally, as a by-product of the analysis of Sec. IIVI we 
obtain scaling results for the interaction of a high-valent 
counterion with a charged plane immersed in a weakly 
correlated ionic liquid. All the results of Sec. IIVI and in 
particular the scaling form (|3TJl . valid for z > \/S, can 
be regarded as exact in such a system. 

Our approach can be easily generalized to more com- 
plicated geometries than the planar one, although the 
practical solution of the PB equation with a test charge 
may be more difficult in these cases. Other natural gener- 
alizations are to consider non-uniformly charged surfaces 
and charged objects in contact with a salt solution. Be- 
yond the TCMF approximation of Eqs. (J2DJ, and 
H2ip . the exact equation 11221) always applies and can be 
a very useful tool for assessing correlation effects near 
charged macromolecules of various geometries. 

We conclude by noting that important questions re- 
main open regarding the infinite planar double layer, 
which is the most simple model of a charged macroion 
in solution. One such issue, on which the present work 
sheds light, is the crossover from a strongly coupled liquid 
close to the charged plate to a weakly correlated liquid 
further away. In particular, the precise functional depen- 
dence of the effective Gouy-Chapman length b on 3 is 
still not known. More simulation results, in particular 
at large distances from the charged plate, and a direct 
evaluation of f(z) from simulation, may be useful in or- 
der to gain further understanding and to test some of the 
ideas presented in this work. Another important issue, 
which has not been addressed at all in the present work, 
is the possible emergence of a crystalline long range or- 
der parallel to the plane at sufficiently large values of 
the coupling parameter. Although plausible arguments 
have been presented for the occurrence of such a phase 
transition at 3 > 3 x 10 4 its existence has not been 
proved. 
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APPENDIX A: MEAN FIELD FREE ENERGY 



APPENDIX B: DERIVATION OF IDENTITY <l23l) 



In this appendix we show how the mean field free 
energy l|18|) is obtained as an approximation to F(zq), 
Eq. (|14(l . We start from a general expression for the 
grand canonical potential of an ionic solution interacting 
with an external and fixed charge distribution er(r). In 
the mean field approximation [20l l27j , 



d 3 r 



1 



8nl B q 2 
- Ae(r)e-^ r >} 



a(r)4>(r) 



(Al) 



where q is the valency of the ions, A is the fugacity, 0(r) 
is equal to 1 in the region accessible to the ions and to 
zero elsewhere (equal in our case to 0(z), the Heaviside 
function), and £1 is given in units of fc^T. Requiring an 
extremum with respect to if, the reduced electrostatic 
potential, yields the PB equation which determines the 
electrostatic potential and the actual value of fl. We 
use equation (|A1I) . which is given in the grand-canonical 
ensemble, because it is widely used in the literature [2(], 
|27| . In Ref. [2(j Eq. (|A1|) is derived systematically as 
the zero-th order term in a loop expansion of the exact 
partition function. 

Inspection of Eq. I|14|) shows that it describes the free 
energy of an ionic solution interacting with an external 
charge distribution having the following parameters, 



1 



valency q 
Ext. potential er(r) = — j= 
Bjcrrum length lg =1 



Z7T 



(A2) 



In the second line (external potential) the first term 
comes from the uniformly charged plate and the second 
term comes from the fixed test charge. Plugging these 
values in Eq. I|A1|I yields, 



= - J d 3 r j-^-(V^) 2 - X9(z)e- 



--^-S(z) + E5(r - z z) 



(A3) 



In order to obtain Eq. (|18fl two modifications are re- 
quired. First, we need to return to the canonical en- 
semble by adding /iTV to f2, where N is the total number 
of ions. Noting that /i = logA and that from charge neu- 
trality qN — — / d 3 rer(r), this modification yields the 
extra term that is proportional to log A in Eq. (|18fl . Sec- 
ond, we note that f2 includes the Coulomb self-energy of 
the charged plane and of the test charge. This infinite 
term does not depend on z and should be subtracted 
from f2 since it is not included in the definition of F(zq), 
Eq. 1H|). 

We finally note that the results of this Appendix can 
also be obtained directly from the canonical partition 
function, as expressed by Eq. I|14|) . 



We would like to evaluate the variation 5Fp^(zo)/5zo, 
where Fpe is given by Eq. I|18|) . Note that ip itself de- 
pends on zq. However the first derivative of Fpb with 
respect to <p(r) is zero. Hence the only contribution to 
SFpb/5zq comes from the explicit dependence on zq: 



SF pb [zq] 
8z 



If d 

- / d 3 r(<^-logA)^— 5(r-z z) 

£2 J OZ 



dp 
~dz~ 



(Bl) 



Z Z 



It is also instructive to derive this identity within the 
exact theory. Equation (|14|l can be written as follows: 



1 f N ~ X 

cxp[-F(z )] = (jv_ !); J II d a riexp(-fl^{ri}) 



where the N-ih charge is fixed at r = zqz: 
AT-l jv-i „ 



(B2) 



Z Z\ 



E 

j>i 



Differentiating with respect to zq yields: 

/ JV-l 



L 3\ 

(B3) 



SF(z ) 
Sz 



dH Z0 
8zq 



- + £ 



^(z Q - Zj) 

i |ri-z z| 3 



(B4) 

where the averaging is performed over all configurations 
of the N — 1 ions with the weight exp(— H ZQ {ri}). This 
quantity is the mean electrostatic field acting on a test 
charge at zqZ. 



APPENDIX C: NUMERICAL SCHEME 

Numerically solving a non-linear PDE such as i|19f) re- 
quires careful examination of the solution behavior. The 
purpose of this section is to explain the numerical scheme 
used in this work, and in particular the parameters re- 
quired to obtain a reliable solution. 



Finite cell 

We solve Eq. <|19|l as a two dimensional problem in 
the coordinates r and z, making use of the symmetry of 
rotation around the z axis. The problem is defined within 
a finite cell of cylindrical shape, shown schematically in 
Fig. 11. The negatively charged plate is at z = and 
ions are only allowed in the region z > 0, marked as the 
gray-shaded region in the plot. 

We impose a boundary condition of zero electrostatic 
field, 



W(p ■ n = 0, 



(CI) 
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R 



where ip n represents the n-th iteration. 




-Z- 



mm 







^max 



FIG. 11: Schematic illustration of the setup used to solve 
numerically the PDE 1191 . The equation is solved in a finite 
cylindrical cell, extending from — z m ia to z ma x in the z axis 
and from to R in the r axis. The charged plane is at z = 
and ions are only present for z > 0. Neumann boundary 
conditions (Vip ■ n = 0) are imposed at the cell boundaries. 
The test charge is at r = zqz 



at the cell boundaries: z = —z m - m , z = z max , and r = R. 
The cell size, as determined by these boundaries, must 
be sufficiently large, as will be further discussed below. 

In the numerical solution it is necessary to solve ip for 
the electrostatic potential at positive as well as negative 
z j3||. Note that a boundary condition such as (jClfl at 
z = would correspond to zero dielectric constant at 
z < 0, while we are interested in the case of continuous 
dielectric constant across the plate. 



Differential equation 

The source term in Eq. (|19fl diverges at z = and at 
r = zqz. We avoid this difficulty by shifting the potential: 



<p = ip + \z\ + 



z z\ 



(C2) 



and solving for tp, which is the potential due only to the 
mobile ions. The equation for ip, 



V 2 ?/> = -47rA6»(z)exp ( —ip - z 



z z\ 



(C3) 



is solved with a Neumann boundary condition for tp, de- 
rived from Eqs. ijClfl and (|C2|) . Note that, unlike <p, ip 
is well behaved at z z. The nonlinear equation l|C3|) can 
be solved by iterative solution of a linear equation (see, 
for example, IHIH ), 



= — 47rA0(z)exp ^—ip n _i 
x [1 - (ip n - V'n-l)] 



z Q z\ 



Grid and solution method 

In the coordinates r,z the cylindrical cell is a rectan- 
gular domain, 

[0,-/?] X [ ^min ? ^max] ■ 

We use bi-cubic Hermite collocation |3(j in this domain 
in order to translate the PDE into a set of linear al- 
gebraic equations on a grid. These equations are then 
solved using Gauss elimination with scaled partial piv- 
oting. Storing the band matrix representing the linear 
equations requires approximately 48 x iV 2 x N z memory 
cells, where N r and N z are the number of grid points 
in the r and z directions, respectively [3(|. Because this 
number can be very large, it is essential to use a vari- 
ably spaced grid in both of the coordinates. We use the 
following scheme: 

r coordinate: In the absence of a test charge, an 
arbitrarily coarse grid can be used in the r direction, 
due to the translational symmetry parallel to the charged 
plane. In our case (where a test charge is present) the 
grid spacing is determined by the distance from the test 
charge, as follows, 



dn 
dr 



r grid 



(C5) 



where n r and r gr ;d are two fixed parameters, while n 
stands for the grid point index and dn/dr is the number 
of grid points per unit increment of the radial coordinate. 
This spacing is approximately uniform up to the thresh- 
old r gr id, whereas for r ^> r gr id it is proportional to 1/r. 
The grid points are then of the form 



'"grid X 







- 1 


X 


exp ( — ] 




\n r J 





(C6) 



In practice r gr id is chosen approximately proportional to 
VS, in order to allow the structure of the correlation hole 
to be represented faithfully. 

z coordinate: In this coordinate the grid spacing is 
influenced by the distance from the charged plate as well 
as the distance from the test charge. We describe sep- 
arately the spacing determined from each of these two 
criteria; the actual grid is obtained by using the smaller 
of the two spacings at each point. 

(i) Distance from the plate: we use a grid spacing pro- 
portional to the derivative of ippb(z): 



1~ r 

dz z+l 



(C7) 



Ignoring, for the moment, the distance from the test 
charge, Eq. i|C7D leads to grid points of the form 



(C4) 



exp(n ■ D) — 1 



(C8) 
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TABLE II: Parameters used in numerical solution of the PDE. 



and 









R 


D 




n z 


fgrid 


n r 


10 4 


10 4 


10 s 


10 5 


0.2 


5000 


5 


100 


5 


10 3 


10 4 


10 s 


10 5 


0.2 


500 


4 


33 


4 


10- 


10 3 


10 4 


10 4 


0.2 


50 


4 


10 


4 


10 


10 3 


10 4 


10 4 


0.2 


5 


4 


3 


4 


1 


4 x 10 2 


4 x 10 3 


4 x 10 3 


0.2 


0.5 


4 


1 


4 


0.1 


80 


800 


800 


0.2 


0.05 


4 


0.2 


4 



where the parameter D is the grid spacing close to the 
charged plate. A similar scheme is used in the negative 
z half-space. 

(ii) Distance from the test charge: we use a form sim- 
ilar to l|C5jl . 



dn 



dz \z- Z \ +2 g rid' 



(C9) 



In practice, the threshold z gr id is chosen proportional to 
H, in contrast to r gr id which is chosen proportional to 



dF(z) 



da 



dF(z) 



dz 



We now use Eq. (1221) to obtain, 



n(0) 



d 



(D4) 



(D5) 



ZqZ 



The second term in this equation is the average electro- 
static force acting on the ions per unit area. This force 
can be separated into two contributions. The first one, 
exerted by the charged plane, is equal to — J dz 1 n(z') — 
— 1 because the plane applies an electrostatic force which 
does not depend on z' and is equal to unity in our rescaled 
coordinates. The remaining contribution to the force, 
exerted by the ions on themselves, must be zero due to 
Newton's third law, leading to the result n(0) = 1. 

The discussion up to now was exact. It also applies to 
PB theory, where the self consistency of the mean field 
approximation ensures that Newton's third law is obeyed. 
On the other hand, within our approximation the force 
exerted by the ions on themselves, 



Parameters 

The parameters that were used to obtain the numerical 
results presented in this work are summarized in Table 2. 
We compared our results with those obtained with (a) 
Increasing z m - ln , z max , and R by a factor of 10; and (b) 
decreasing the grid spacing by a factor of 2, both in the r 
and in the z coordinates. The influence of these changes 
was found to be negligible on all the data presented in 
this work. 



APPENDIX D: CONTACT THEOREM 

In this appendix we derive the contact theorem [23j in 
a way that highlights the reason why it is not obeyed in 
our approximation. We start from an exact expression 
for the free energy, 



/>oc 

F = -log / dz' exp[-F(z')} 

J a 



(Dl) 



where the charged plate is at z — a. This plate position 
can be chosen arbitrarily, hence dF/da = 0: 







OF 
da 



n(0) 



' dz'n(z')?f(z>) (D2) 
dz 



n(z) [f{z) - 1] 



is not zero. This inconsistency can be traced to a more 
fundamental inconsistency which is briefly described be- 
low. 

The probability to find the test charge at r = zz and 
a mobile ion at r' is proportional, in our approximation, 
to 

n(z)n(z')g(r, r') = n(z)exp [— </?(r'; z)] (D6) 

In the exact theory the probability to find two ions at 
r and r' must be symmetric with respect to exchange 
of r and r'. On the other hand the correlation function 
g(r, r'), as defined above, is not symmetric. In other 
words, the ion-ion correlation function in the TCMF 
model is not symmetric. 



APPENDIX E: SMALL S EXPANSION 

The recovery of mean field results at small 2 was 
demonstrated and explained in Sec. [n] Here we derive 
this result formally as an expansion in powers of S. The 
advantage of this formal expansion is that it allows us 
to find also the first order correction to the PB profile 
within our model. 

We expand SFpb, Eq. 1)18(1 . up to second order in 3: 



where we used the relations 



ZF PB {z ) =F a + SFi(zo) + S 2 F 2 (z ) + 



(El) 



The zero-th order term, Fq 7 does not depend on zq and is 
the PB free energy of a charged plane in contact with its 



17 



counterions, without a test charge. In order to evaluate 
the following terms, we also expand ip in powers of S: 

ip(r; z Q ) = tp Q (v) + 5<^i(r; z ) + E 2 tp 2 (r; z ) + ■ ■ ■ (E2) 

To zero-th order we have from Eqs. I|18|l and i|19[l : 



F Q = |d 3 r|-i-(V^ o ) 2 -A0(z)e-^( z )| 



where 



dz 2 



-47rA6»(z)e" 



(E3) 



(E4) 



is the potential due to counterions in the PB approxi- 
mation. The first order term in the free energy, F±, is 
found by expanding equation l|18(l in S. This expansion 
includes two contributions, the first from tpi and the sec- 
ond from the explicit dependence on S in Eq. (|18fl . The 
first contribution vanishes because tpo is an extremum 
of the zero-th order free energy, leaving only the second 
contribution: 

Fi(zo) = J d 3 rtpo(r)5(r - z i) = tpo(z ) (E5) 

Returning to our approximation for n(z), given by 
Eq. JU, we find that: 



n(z) = iexp[-F PB (z)] = 4 ex P 



Z 
1 



-— - ipo(z) 



exp [-ipo(z)} 



(E6) 



where Zo is found from the normalization condition (|21|l . 
To leading order in 5, n(z) is equal to the PB density 
profile, as expected: 

n(z) = n PB {z) = exp[-tp (z)] = ; - : ^ 2 (E7) 
Z {z + iy 

where Zq is a normalization constant. The next order 
term in the expansion of / can be found on similar lines 
as F\ (z) , and is equal to 



^2(20) = -S(pi(z z; z ) 



(E8) 



where Sip is the difference between the first order cor- 
rection to tp and the bare potential of the test charge: 



The function Stpi (r) arises also in the systematic loop ex- 
pansion of the free energy around the mean field solution 
|20|. Its value at r = z z is given by p0| : 



5tpi(z z; z a ) = g(z ) = 



2(z + l) 2 
{le^-^E, l(l-i)z }(l + iz f 

-ie^ Z0 E 1 [(1 + i)zo] (1 - iz f - 4z } (Ell) 

where E\ [x] is the exponential-integral function |25| . Us- 
ing Eqs. (|20|l and (|E8|I we find that up to first order in 
S the density profile is given by: 



n{z) = ^pb(^) + Sni(z) 



where 



ni(z) 



n PB (z) 



(E12) 



(E13) 



In this expression g(z) is given by Eq. (|E 1 lfl and iVi is 
obtained from the normalization condition 121(1 : 



Ni = 



dz n PB (z)g(z) 



-0.3104 



(E14) 



Note that this is different from the exact expression for 
the first order correction in S [37]], which is obtained in 
the loop expansion and is not reproduced here, but is 
shown in Fig. 12. 

Figure 12 shows n\{z) as defined by Eq. IE 13(1 (solid 
line). The symbols show the correction to ripe calculated 
numerically from TCMF for 3 = 0.1 and scaled by 1/H = 
10. At this small value of S the linearization provides a 
very good approximation for the correction to n PB {z). 

The dashed line shows the exact first order correction 
in 5 to the ion density, obtained from the loop expan- 
sion. Comparison of the solid and dashed lines shows 
that the TCMF model does not capture correctly the ex- 
act first order correction. In particular, ni(0) is different 
from zero in our approximation; in the exact correction 
ni(0) = as it must be due to the contact theorem. It 
is important to realize that although the exact first or- 
der correction is useful for values of S of order unity and 
smaller, the TCMF has a much wider range of validity 
for 5 > 1. 



Proof of Equation \3J$ 



5yi{r) = <pi(r) - 



r - z z 



(E9) 



The first order term in the expansion of tp, tpi(r; zq) is 
the solution of the differential equation: 

[V 2 - 47rAe- V0 ] tp 1 = -4nS(r - z a z) (E10) 



Our purpose here is to prove the first equality of 
Eq. El, 



h(z) 



1 dgjz) 

2 dz 



(E15) 



where the electrostatic field acting on a test charge is 
— (/pb(z) + 2/i (z) + ...), i.e., fi(z) is the first order 
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0.14r 



N 




3 
Z 



FIG. 12: First order (linear) correction in S to n(z), as 
obtained from the test charge mean field approximation, 
Eq. l|E13fl (solid line), compared to the exact first order cor- 
rection calculated using a loop expansion (20ll (dashed line). 
The symbols show [n(z) — npB(z)]/H calculated numerically 
in the test charge mean field approximation with S = 0.1. 



term in 3. In order to do this, let us consider the correc- 
tion to the mean field potential due to an infinitesimal 
point charge of magnitude H that is placed at r = zz. 
We designate this correction, evaluated at the point r', 
as G(r,r'). This Green's function is found by solving 
Eq. i|E10() which reads, with a slight change of notation: 



G(r, r') = -47r£(r - r') (E16) 



V 2 , - 4 7 rAe- ¥ ' o(r ' ) 



The electrostatic field acting on the test charge is then 
-/pb(z) - 2/i(z), where 



d 



h{z)=—G{vy) 



= IA G (r r )- ldg(z) 
r ' = r 2 8z 1 ' ' 2 dz 



(E17) 

and g(z) is defined in Eq. IjElip . In the second step we 
used the symmetry of G(r,r') to exchange of r and r', 



which follows from the fact that the operator acting on 
G(r, r') in Eq. I|E16|I . as well as the right hand side of 
that equation, are symmetric with respect to exchange 
of r and r'. 



APPENDIX F: MEAN FIELD EQUATION AT 
LARGE z 



We start from the exact identity (|43() and would like 
to evaluate f(z) for a test particle placed at sufficiently 
large z, assuming also that 3 is large. The mean field 
electrostatic force acting on the particle is given by 

f-Z />oo 

/mf(z) = 1-/ dz'n{z') + J dz'n(z') (Fl) 

where the first term on the right hand side is the con- 
tribution of the charged plane, the second term is the 
contribution of ions between the plane and the test par- 
ticle, and the third term is the contribution of the other 
ions. Eq. (|F1I) would describe the exact force acting on 
the test particle had it not had any effect on the distri- 
bution of the other ions in the system. We need to add 
to this force the contribution due to the influence of the 
test charge on the other ions. 

Due to the exponential decay close to the plate the ion 
layer further than z = \/3 is very dilute. Hence it makes 
sense to include in the correlation-induced force only a 
contribution from the ions close to the plate. Estimating 
this contribution as aB./z 2 we conclude that 



dlogn(z) 
dz 



-f( z ) = ~/mf(z) 



(F2) 



Differentiation of this equation with respect to z yields 
Eq. gSJ): 



d 2 logn(z) 2aZ 
dz 2 - Zn y z l+ Z 3 



(F3) 



[1] Electrostatic Effects in Soft Matter and Biophysics, 
edited by C. Holm, P. Kekicheff, and R. Podgornik 
(Kluwer Academic Publishers, Dordrecht, 2001). 

[2] W.M. Gelbart, R.F. Bruinsma, PA. Pincus, and V.A. 
Parsegian, Phys. Today 53, 38 (2000). 

[3] R. Podgornik and B. Zeks, J. Chem. Soc. Faraday Trans. 
2 84, 611 (1988). 

[4] M.J. Stevens and M. O. Robbins, Europhys. Lett. 12, 81 
(1990). 

[5] R. Kjellander, T. Ak esson, B. Jonsson, and S. Marfielja, 

J. Chem. Phys., 97, 1424 (1992). 
[6] B.Y. Ha and A.J. Liu, Phys. Rev. Lett. 79, 1289 (1997). 
[7] R.R. Netz and H. Orland, Europhys. Lett. 45, 726 (1999). 
[8] D. B. Lukatsky and S. A. Safran, Phys. Rev. E. 60, 5848 



(1999). 

[9] V.I. Perel and B.I. Shklovskii, Physica A 274, 446 (1999); 

B.I. Shklovskii, Phys. Rev. E. 60, 5802 (1999). 
[10] R.R. Netz, Eur. Phys. J. E. 5, 557 (2001). 
[11] L. Guldbrand, B. Jonsson, H. Wennerstrom, and P. 

Linse, J. Chem. Phys. 80, 2221 (1984). 
[12] N. Gr0nbech- Jensen, R.J. Mashl, R.F. Bruinsma, and 

W.M. Gelbart, Phys. Rev. Lett. 78, 2477. 
[13] A.G. Moreira and R.R. Netz, Eur. Phys. J. E. 8, 33 

(2002) ; A.G. Moreira and R.R. Netz, Europhys. Lett. 
52, 705 (2000). 

[14] M. Deserno, A. Arnold, and C. Holm, Macromol. 36, 249 

(2003) . 

[15] R. Kjellander, S. Marcelja, and J. P. Quirk, J. Colloid Int. 



19 



Sci. 126, 194 (1988). 
[16] V.A. Bloomfield, Biopolymers 31, 1471 (1991). 
[17] D.C. Rau and A. Parsegian, Biophys. J. 61, 246 (1992); 

D.C. Rau and A. Parsegian, Biophys. J. 61, 260 (1992). 
[18] P. Kekicheff, S. Marcelja, T.J. Senden, and V.E. Shubin, 

J. Chem. Phys. 99, 6098 (1993). 
[19] E. Raspaud, M.O. de la Cruz, J.L. Sikorav, and F. 

Livolant, Biophys. J. 74, 381 (1998). 
[20] R.R. Netz and H. Orland, Eur. Phys. J. E. 1, 203 (2000). 
[21] R. Kjellander and S. Marcelja, J. Chem. Phys. 82, 2122 

(1985); R. Kjellander, J. Chem. Phys. 88, 7129 (1988); 

R. Kjellander and S. Marcelja, J. Chem. Phys. 88, 7138 

(1988). 

[22] For a review, see D. Andelman, in Handbook of Physics of 
Biological Systems, edited by R. Lipowsky and E. Sack- 
mann (Elsevier Science, Amsterdam, 1994), Vol. I, Chap. 
12, p. 603. 

[23] S.L. Carnie and D.Y.C. Chan, J. Chem. Phys. 74, 1293 
(1981). 

[24] A.G. Moreira and R.R. Netz, private communication. 

[25] M. Abramowitz and LA. Stegun, Handbook of Mathemat- 
ical Functions (Dover Publications, New York, 1972). 

[26] L. Blum and D. Henderson, in Fundamentals of Inhomo- 
geneous Flmds, edited by D. Henderson (Marcel Dekker, 
New York, 1992), Chap. 6, pp. 239-276, and references 
therein. 

[27] Y. Burak and D. Andelman, Phys. Rev. E. 62, 5296 
(2000). 

[28] D.L. Harries, S. May, W.M. Gelbart, and A. Ben-Shaul, 

Biophys. J. 75, 159 (1998). 
[29] S.L. Carnie, D.Y.C. Chan, and J. Stankovich, J. Colloid 

Interface Sci. 165, 116 (1994). 



[30] E.N. Houstis, W.F. Mitchell, and J.R. Rice, ACM Trans. 
Math. Soft. 11, 379 (1985). 

[31] A related dimensionless parameter is L, conventionally 
defined for the two dimensional one component plasma 
(see Ref. 0). This parameter is related to E as follows: 

s = 2r 2 . 

[32] The numerical differentiation of log[n(z)] leads to large 
error bars because n(z), as obtained from the simulation, 
is noisy. In principle f(z) could be evaluated more accu- 
rately during the simulation run by direct use of Eq. (1241 , 
i.e. by averaging the electrostatic force acting on ions as 
function of their distance from the plane. 

[33] It is also possible to find n(z) using Eqs. (Unj-HD, but 
integrating Eq. 1251 is numerically more accurate. 

[34] In principle, a more accurate evaluation of the pair distri- 
bution function, possibly capturing its oscillatory nature, 
may be obtained from an approach similar to the TCMF 
with two test charges instead of one. 

[35] In order to improve over TCMF results it may by use- 
ful to treat mobile ions close to the plane as a two di- 
mensional layer, while going beyond the mean field ap- 
proximation of the TCMF model in their treatment. The 
response to the test charge at r = zoi, may then help un- 
derstand the decay of n(z) in the full three dimensional 
problem, through Eq. 1251 . Such a calculation is beyond 
the scope of the present work. 

[36] This situation is different from that of the PB equation 
with no test charge, in which the electrostatic potential 
depends only on z. The boundary condition at z = 0, 
for the case of zero dielectric constant at z < 0, is then 
sufficient in order to solve for the potential at z > 0. 

[37] See Eq. (64) and Fig. 4 in Ref. Hj. 



